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This report documents the development of a two dimensional finite element based 
numerical model for efficient characterization of the Hall thruster plasma dynamics in the 
framework of multi-fluid model. Effect of the ionization and the recombination has been 
included in the present model. Based on the experimental data, a third order polynomial 
in electron temperature is used to calculate the ionization rate. The neutral dynamics is 
included only through the neutral continuity equation in the presence of a uniform neutral 
flow. The electrons are modeled as magnetized and hot, whereas ions are assumed 
unmagnetized and cold. The dynamics of Hall thruster is also investigated in the presence 
of plasma-wall interaction. The plasma-wall interaction is a function of wall potential, 
which in turn is determined by the secondary electron emission and sputtering yield. The 
effect of secondary electron emission and sputter yield has been considered 
simultaneously. Simulation results are interpreted in the light of experimental 
observations and available numerical solutions in the literature. 

Nomenclature: 

B, B Magnetic field, gauss 

E Electron charge, coulomb 

Ei Ionization potential, eV 

E, E Electric field, V/m 

j, J Current, mA 

L Differential operator 

M Mass, kg 

M Mass matrix 

N Number density, m ' 3 

N Basis function 

R Radial direction 

R Solution residual 

S Assembly operator 

T Time, s 

T Temperature, eV 

U, U State variable 

V, V Velocity, m/s 
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f 


w 

Weight function 

Xe 

Xenon 

Z 

Axial direction 

Z 

Ionicity 

Greek 

V 

Collision frequency, s' 1 

</> 

Potential, volt 

r 

Flux of the propellant, m'V 1 

Q 

Solution domain 

a 

Ionization cross-section, m 2 

3 

Implicitness 

A 

Step 

Subscripts 

*,0 

Reference value 

c 

Charge exchange 

D 

Discharge 

e 

Electron 

i 

Ion 

k 

Degree of interpolation polynomial 

H 

Hall current 

n 

Neutrals 

r 

Radial component 

t 

Thermal velocity 

2 

Axial component 

a 

Electron or ion 

0 

Azimuthal component 

T 

Time stepping index 


Superscript 


s 


h 

Discretization 

P 

Iteration index 

0 

Neutrals 

+ 

Singly ionized 

++ 

Doubly ionized 


I. Introduction 

Hall effect thruster (HET) is a highly effective on-board propulsion device popular in many LEO and 
GEO satellites. 1 It is also called “closed drift” thruster (CDT) due to the azimuthal drift of electrons. 
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This is common to variants of such thrusters, e.g. stationary plasma thruster (SPT), thruster with anode 
layer (TAL) etc. The SPT thruster is a coaxial device that consists of four main parts: the anode, which 
serves as a propellant distributor; an annular acceleration channel made of boron nitride; a magnetic 
unit; and a hollow cathode (Fig. 1). The plasma column is contained within two coaxial dielectric 
cylinders that constitute the discharge channel, with the anode at one end of the channel and the exit at 
the other end of the channel. The discharge is created between the anode of the thruster and an external 
hollow cathode located downstream of the channel exit. The magnetic system consists of a series of 
electromagnetic coils employed inside the inner cylinder and outside the outer cylinder, and pre- 
dominantly radial field with a maximum just upstream of the channel exit. The electrons from the 
cathode enter the chamber 
and are subject to 
azimuthal E*B drift. The 
electrons in the closed 
drift undergo ionizing 
collisions with the 
propellant gas. While the 
magnetic field is strong 
enough to capture 
electrons in an azimuthal 
drift, it is not strong 
enough to contain the 
resulting ions, which are 
essentially accelerated by 
the imposed axial electric field. The suppression of axial electron mobility by the imposed radial 
magnetic field, while leaving ion mobility unaffected, enables the plasma to support an electric field 
with a potential difference close to the applied voltage. The ions are accelerated to kinetic energies 
within 80% of the applied discharge voltage. 2 

Present day Hall thrusters offer specific impulses over 2400 s, thrust over 1 N, and power exceeding 50 
kW at efficiencies close to 60%. 3 However, the commercial exploitation of Hall thrusters imposes a 
stringent constraint of trouble free operation for more than 8000 hrs. 4 The physics inside the Hall 
thruster has to be reasonably well understood in order to make any significant change in efficiency 
without compromising the lifetime. This is a challenge, as the choice of thruster size requires an 
optimum selection between efficiency and lifetime. 5 Despite significant numerical and theoretical 
advances of the recent past, we lack an adequate numerical model to describe critical regions of a Hall 
thruster plasma dynamics in a self-consistent fashion. 6 ' 7 This effort was focused towards developing a 
consistent formulation for practical SPT simulation. 

Numerical simulation of the plasma dynamics of a Hall thruster has been carried out recently by several 
authors in the framework of the hybrid as well as the fluid models. 8 ' 28 The one-dimensional (ID) fluid 
model of the partially ionized plasma incorporating the neutral dynamics and the effect of the plasma- 
wall interaction has been documented recently. 26 ' 28 This report presents the numerical and theoretical 
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development a 2D, three-fluid, partially ionized plasma model in order to investigate the effect of 
ionization and recombination on the dynamics of the Hall thruster. The neutral dynamics is included in 
the present work since without neutral dynamics the effect of ionization and recombination cannot be 
studied satisfactorily. The self-consistent 2D, three-fluid formulation of the bounded thruster plasma is 
the novel feature of this report. To the best of our knowledge, such a simulation has not been reported in 
the literature. 

Numerical novelty includes the utilization of sub-grid embedded (SGM) finite elements, 29 ' 30 for 
convergence and stability of the solution. It is based on a non-linear, non-hierarchical, high-degree 
Lagrange finite element basis for use in a discretized approximation. SGM elements utilize local mesh, 
velocity and diffusion parameters to modify the dissipative flux vector (second derivative) terms in the 
equation. For the hyperbolic equation, a second derivative artificial diffusion term with a vanishing 
coefficient is added. The theory employs element-level static condensation and eigenvalue analysis for 
efficiency, nodal-rank homogeneity and essentially non-oscillatory solution. Unlike traditional upwind 
methods, however, non-linear SGM does not introduce any unnecessary diffusion to distort the solution. 


The numerical model and the simulation results are presented in the subsequent sections. In section n, 
we discuss pertinent theoretical issues. In section III, the solution algorithm is described. The numerical 
results are documented in section IV. Finally, section V contains conclusions and future work. 

II. Theoretical Issues 

The dynamics of a partially ionized, thruster plasma is quite complicated, 1-4 ’ 8 ' 13 ’ 16 ’ 18 ' 21 ’ 27 ’ 28 as several 
elastic and inelastic processes may occur simultaneously. However, not all processes are equally 
probable. For example, momentum exchange between electron-electron and ion-ion will not be 
important in comparison with the electron-ion momentum exchange as the relative drift between similar 
particles is small in comparison with the drift between electrons and ions. The collisions between 
electron-neutral, electron-ion, and ion-neutral play an important role. The plasma-neutral collision 
usually determines the kinetics of the motion. 


(a) Ionization: Electrons collision with the neutral atom is the main source of ion production in 
propulsion plasma. The rate of ion production in plasma is determined by the total cross section of the 
process. 


Sion* = (^( 7 )) = 


( 6 ) 


where, for process constant k = <Va'(V)> , the averaging is done over the velocities of the electrons 
whose energy is sufficient for ionization mV 2 12 >£,. The ionization source term, which takes into 
account all the above processes, is 


Siangan = 


( 7 ) 


where 0+, 0++ represents the transition from neutral to singly and doubly ionized state respectively 
(Xe°— >Xe + ’ and Xe°— >Xe ++ ) and 1++ represents the transition from singly to doubly ionized state 
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(Xe + ->Xe ++ ). A general third order temperature dependant polynomial can be fitted to the experimental 
value of ionization rate k t = &, 0+ + &, 0++ + k/ + . The matrix form is 


'tit* ' 


l0++ 




J++ 


' 1.9435xl0' 5 
-3.0352xl0" 5 
-2.117xl0' 5 


-0.0068 0.6705 -1.6329 
0.0024 0.0515 -0.1431 
0.0022 - 0.0119 0.0161 


f r T 3 > \ 

L e 

Tt 

t 


xlO 




( 8 ) 


The sum of all three ionization rates can be represented as, 

k, = (-3. 2087x1 O' 5 T]~ 0.00227/ + 0.7 10 \T e - 1.76) xl O' 14 (9). 

Above estimate of ionization rate is based on the Maxwellian distribution function. 


fbl Recombination: The plasma-wall interaction leads to the recombination of the plasma particles at the 
wall. Furthermore, the recombination in the presence of a neutral body (wall) is important at the low 
degree of ionization. The recombination coefficient a can be approximated as 31 
a = 1.09x1 0" 20 n e T* n m 3 /s. Then the recombination rate can be written as 

S r ecom = -1-09x1 0 _2 ° T' 9,2 n 3 m 3 /s , where assuming quasi-neutrality «,■ has been replaced by n e . 


(ct Charge-exchange process: Charge exchange is related with the transfer of one or more electrons 
between an atom and an ion. Slow propellant ions are created due to resonant charge-exchange 
collisions between the fast “beam” (current) ions and slow thermal neutrals. The spatial volumetric 

production rate is given by S CEX =n ri n i (v l c(y i )}, where relative collision velocity is taken to be the ion 


■Xe^+XeL, XC+XcJL 


>Xe~ +XeL, Xe~ +Xe» 0W 


'Xe^+Xel 


velocity. 

Xe^rt +Xe°| 0W -T/w i|ow ^w fast T^vw jl0W — r '‘•''slow 1 '■"'fiist’ -^fut ^ ^*- w 8low ^ ^■‘ v 8 low 

The process can be important for creating slow ions. The collisional cross-section for the above 
processes are comparable 31 cr(Xe-Xe + ) ~ 4.38x1 0" 19 m 2 and cr(Xe-Xe ++ ) ~ 4.98x1 0" 19 m 2 . Thus, it 
would appear that all the above processes are equally important. However, in the last process stripping 
of an electron may require energy exceeding 1 keV 33 and thus, the last process is ignored. The cross 
section for Xe-Xe for example is given by 34 


o-(Xe - Xe + ) = [a - b log I0 (Am )] (E, jE H ) 


■ 1.5 


xlO" 20 m 2 . 


( 9 ) 


where a = 181, b = 21.2, £, = 12.13 eV - xenon ionization potential and Eh - 13.6 eV - hydrogen 
ionization potential. For a relative velocity Am between 10 and 2x10 m/s, the charge exchange cross- 
section varies between 10' 20 to 10' 19 m 2 . 


Having delineated some of the important physical processes in the partially ionized thruster plasma, we 
shall now give the basic set of equations that describes the dynamics of the process under investigation. 
We shall assume that the ions are unmagnetized, since for typical parameters of a thruster plasma viz., 
magnetic field B ~ 200G and ion velocity 4xl0 3 m/s the gyration radius of ions are about 0.1 m, which 
is much larger than the size of the thruster (0.02 - 0.03 m). In the present simulation, the maximum 
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value of magnetic field, near the exit, is 242 G. Therefore, the effect of magnetic field on the ion 
transport will be ignored. Further, the pressure term in the ion momentum equation can be ignored as the 
thermal energy of the ions is much smaller than their kinetic energy i.e. Ti«mjVj 2 . Note that owing to 
the small inertia, electron response time is much faster than the ion response time. As a result, electrons 
will attain the steady state faster than the ions. Keeping this in mind, electron momentum and energy 
equations are solved at steady state, whereas for ions and neutrals, a set of time-dependent continuity 
and momentum equations are simultaneously solved. 


An annular cylinder can adequately characterize Hall thruster geometry. Ignoring any variation in the 
azimuthal (0) direction, we shall take a two-dimensional axisymmetric (r, z) representation of the 
thruster. The following set of equations written in the component form is used to describe the dynamics 
of thruster plasma with the azimuthal ion velocity, and dependence on the azimuthal angle, d/dO, set 
to zero. The continuity equation for the electron and the ions are, 


dn„ Id. . d(n„V„) 

— + (r»X r )+- — = S, 

dt rdr a ar dz 


( 6 ) 


where S = S recom b+ Si 0 ni 2 +Scex ■ Assuming that the neutral velocity has axial component only, the neutral 
continuity equation is 

~ - kfn,n, -C*",".- (7) 

dtdz 

The right hand side of the continuity equations (6)-(7) represents source and sink term due to ionization 
and recombination. The sink term of the neutral continuity equation (7) is slightly different than the 
ionization source term of die plasma continuity equation (6) due to the presence of the processes like 
Xe + — >• Xe 4 ^ in the S ioniz . 

The ion momentum equations are, 
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The factor 0.5 before ion-neutral collision term comes from reduced mass m, m„ /(m,+w„)am,/2. 
The electron momentum equations are, 


K^+K%-£=—^—($+YA)^-K)-y„<K-rd- 

cr a r nyi e ar v\ [n^ 

-Vo-Oi+O^. 


K. 


yE, +v Ei V - V - 

* dr a dz r 


4 | Vrf _ 


\ n eJ 


( 10 ) 

( 11 ) 
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The electron energy equation is, 

l(^ +v -' vr 0 = " r ‘ v ’ v ‘“^' + ^ +K ' +2 ^' (r " +r - )+ ’'“^ 

+*-.(^+^,)+-.(^+^.)+3 i >' J a;-r,)+3iv.(r,-r,)+^|4r,+o£,|. 

m i m n «A 2 ) 

where v = v d + v n . In equations (6)-(13), V e , F are the electron and ion velocities, respectively. E and 


(13) 


B are the electric and magnetic fields, and n ; is the number density of the y-th particle with j = e and i, 
Vei, Ven are the electron-ion and electron-neutral collision frequencies respectively, v c is the ion charge- 
exchange collision frequency, e is the electron charge, Z is the ionicity, and the value of a is between (2 
-3). 19 The electron dynamics are determined by the pressure gradient, electric and magnetic forces, and 
the collisional exchange of momentum in equations (10)-(12). The convective term in these equations 
retains the effect of the electron inertia. 


The contribution of the electron inertia is small and on this basis, its effect on the plasma dynamics is 
generally ignored. However, in the regions of sharp flow gradients, the effect of the convective term 
may become finite, and therefore the convective term is retained in this formulation. Similarly, since 
collision time scales are much slower than the electron-cyclotron gyration time scale, one may ignore 
elastic and inelastic collision terms in comparison with the Lorentz force term VxB in the momentum 
equation. Such an approach will exclude the dynamics of the momentum exchange as well as the effect 
of ionization and recombination, severely limiting the applicability of the model to the thruster plasma. 
Therefore, all the collision terms are retained in the electron momentum equations (10)-(12). 


Equation (13) includes the effect of Joule heating, a contribution due to exchange of random thermal 
energy and due to ionization and recombination. The simulation of Katz et. al 15 shows that the average 
ion energy remains nearly constant in the channel. Therefore, we also assume in our model that the ion 
energy is constant. 


In closure, equations (6)-(13) are supplemented with the simple form of Ohm’s law, 
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Vri'T, (14) 


Equation (14) yields the relationship between the current and strength of the electric field in the plasma. 
The electric field has been written in the moving frame of the plasma, i.e. EsE+VxB. 


Before numerically solving the above set of equations (6)-(14), we normalize the physical variables as 
follows. Temperature T is normalized using the first ionization potential of Xenon, T» = E\ = 12.1 eV. 
All other dependent variables can be normalized from: 
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V.= I— = 2xl0 3 m/s; n. = — = 0.5xl0 17 m' 3 ; v.=c.n.V. = of. □ 2.xl0 5 s' 1 

Vm, V, 


where a.=o 0 /— . The fundamental length scale /o is defined in terms of characteristic velocity and 
V m < 

collisional frequency, /o=V*/v*. The time scale is to = l/v». 


In order to numerically solve the formulation (6)-(14), proper initial and boundary condition 
specifications are necessary to make the problem well posed. The axial velocities of electrons and ions 
are not fixed at the inlet. Under typical conditions, next to the anode, a plasma sheath (typical width ~ 
Debye length) forms and ions must flow into the sheath from the quasi-neutral region. The axial velocity 
is near zero close to anode and then begins to rise at the edge of the acceleration zone, reaching a 
maximum velocity beyond the exit. 22 Such flow behavior has also been observed in the classical nozzle 
problem, where flow changes smoothly from subsonic (in the narrow region) to supersonic flow in the 
divergent region. Therefore, a sonic point, where the flow velocity equals the characteristic speed of the 
medium, is always expected at the exit. We shall impose the “choked exit” 19 boundary condition for the 
ion axial velocity. Furthermore, axial electron flow will be assumed inward at the channel exit. We 
impose zero radial velocity for the electrons and ions at the inlet and leave it floating elsewhere. At the 
inlet, the plasma density is equated to the reference value. The neutral density at the inlet depends upon 
the mass flow rate. In die present calculation, the mass flow rate is ~ 6 mg s' 1 , a value relevant to a 
thruster operating at 1.5 kW power level. 2,14 The corresponding neutral density is rt„ ~ 10 18 tn 3 . The 
homogeneous Neumann condition is imposed on the electrostatic potential at the inlet. Since at the 
cathode, the potential is zero, a vanishing potential is assumed at the outlet. For ion density, a 
homogeneous Neumann condition is assumed. At the upstream boundary (thruster inlet plane), we 
specify an electron temperature of T e = 5 eV, that is close to the experimental data. 6,25 


In a typical Hall thruster experiment, the radial magnetic field is dominant compared to the axial field. 
The radial magnetic field decreases from a typical maximum of about 200 G near the channel exit to a 
much lower value (-30-40 G) near the anode, though higher values of radial magnetic fields (~ 350 G) 
have been utilized recently for P5-class of thrusters 2 The presence of this radial magnetic field inhibits 
electron flow to the anode, and in the process considerably enhances the ionization due to electron 
impact. In the presence of a very strong radial magnetic field and in the absence of any collision, the 
axial electron current may be completely inhibited. Thus, the current is mainly carried by the ions j z ~ j). 
Assuming a quasineutral (n, « n e ) plasma, the Hall current per unit radius is 

J H ~en e dr « en t fa IB (15) 

where fa is the discharge potential. Note that the discharge potential is the sum of the column potential 
drop <f>, the cathode fall potential, and the possible potential drop in the plasma region next to the 
exhaust and outside the cylinders. The corresponding Hall current density can be expressed as 

j H = en e E I B» en e E z / B r (15a) 

For en e Vi, we have 
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(16) 


J H ~ MJBV.m j.^^kr 

H JtTd . 


where the maximum ion velocity is F;=(2e$/w,-) 1/2 . Clearly, for a given magnetic field, J^jr ^<j>a For an 
efficient operating system, current is carried by the ions and ions attain maximum velocity. Thus we 
shall anticipate that the ensuing potential distribution (and the resultant accelerating electric field) will 
be in the region of maximum magnetic field strength. 


Ila. Plasma-Wall Interaction 

The walls of the discharge chamber of a stationary plasma thruster (SPT) are commonly made of 
composite ceramic materials, namely, boron nitride and silicate oxide. Among many reasons limiting 
the efficiency and lifetime of a Hall thruster, the most critical is the wear of the surface layer of the 
ceramic walls. The wall erosion of the thruster occurs due to the plasma-wall interactions. The coaxial 
wall of a thruster develops the non-uniformities due to sputtering, re-deposition, cracking, etc. Further, 
sputtered material may contaminate the spacecraft surface and affect the working parameter 
optimization. Although the lifetime issues are critical to its design, many physical aspects in thruster 
plasma are yet to be understood. The lifetime of an on-board Hall thruster is expected to exceed several 
thousand hours. This complicates the experimental investigation and numerical prediction of the wall 
wear as several parameters come into play during the operational lifetime of the thruster. This results in 
the lack of reliable data on the sputtering yield under operational conditions. 


If the ion impact energy is sufficiently large, the impact ions may cause severe sputtering of the target, 
which is an undesirable side effect. The dielectric walls, anode and hollow cathode walls may develop 
the non-uniformities due to sputtering, re-deposition, cracking, etc. Further, sputtered material may 
contaminate the plasma. This can significantly affect the performance of the HET. The erosion of the 
wall can take place due to the ion bombardment (classical erosion) as well as due to the near wall 
electric fields (anomalous erosion). Whereas, ion bombardment can give rise to small-scale 
prominences mostly across the incident ions, the “anomalous erosion” is manifested as a periodic 
structure oriented along the ion flux with a period of the order of electron Larmor radius 13 , indicating, 
sputtering due to electrons. The difference between the sputtering caused by the atoms and ions 
interaction with the wall is small and thus, it is sufficient to study the ion induced sputtering of the wall. 
This ejection normally occurs when the lattice particle receives sufficient energy from the incoming 
particle to overcome the binding potential of the solid. The minimum bombarding energy needed for 
sputtering is called the threshold energy. Experiments 36 ' 38 indicate that the sputtering yield does not 
depend on the angle of incident ions and it varies linearly with ion energy. Present report employs a ID 
model described below to study the effect of sputtering and secondary emission on the acceleration 
process inside the channel. It is anticipated that the result will provide the basic insight of the underlying 
physics of the plasma-wall interaction. 


Also, the secondary electron emission (SEE) is an important issue. Secondary electron emission is the 
emission of electrons from the wall due to the electron bombardment. A high-energy primary electron 
enters the solid and dissipates its energy. Some of this energy goes into the creation of excited electrons 
and some of these electrons called secondaries escape from the wall. The “intermediate” energy 
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electrons are responsible for the ejection of secondary from the wall. 13 Since the ions and electrons have 
opposite charge, the emitted secondary electrons are accelerated away from target in the same electric 
field that accelerates the ions towards the target. This leads to considerable power loss since part of the 
power goes into the electrons. 38 ' 39 Understanding sheath dynamics is also important in electric 
propulsion devices. Wall effects can significantly alter the dielectric wall characteristics of Hall 
thrusters, while impact of high energy ions may deteriorate the performance of electrodes in MPD 
thrusters. The design of the future nuclear device must deal with the problem of sputtered wall material. 
Sputtered wall material, for example, may be ionized in the scrape-off layer and possibly transported 
into the core plasma, or may be redeposit immediately if the ionization occurs inside the sheath. 

Most of the commonly used plasma confining materials have SEE coefficients near or above 0.9 at 
moderate plasma temperature. For example, boron (15 eV), carbon (12 eV), aluminum (47 eV) etc. have 
SEE which can reach up to O.9. 40 The effect of the cold SEE on the sheath has shown to reduce the 
sheath potential significantly. 8 In fact the onset of space-charge saturation when the electric field is 
reversed near the wall surface, has been found to take place within the plasma sheath when SEE 
coefficient is ~ O.9. 8 Earlier experimental studies of the effect of SEE on plasma sheath potential showed 
that the cold electron emission affect bulk plasma properties. 9,41 

The thruster plasma is partially ionized gas, consisting of electrons, ions and neutral xenon particles. In 
such partially ionized plasma, elastic and inelastic process takes place simultaneously. The elastic 
collision involves only exchange of momentum and energy between colliding particles whereas inelastic 
processes like ionization, recombination, charge-exchange collision, plasma-wall interaction, secondary 
emission, sputtering etc. can be responsible for redistributing the number density of the particles along 
with its momentum and energy. Not all processes are equally probable. We describe the processes 
important for the channel dynamics besides ionization and recombination. 


(a) Electron-ion and plasma-neutral collisions: Because of the long-range nature of the Coulomb force, 
plasma particles can be deflected over the Debye length X D . The electron-ion collision frequency is 


= 4a/2k e 4 n i L e _ 
1 3 Tf 


K 

r \ 

n i 

fO 

3>/2tc 

UJ 

^ n c^Dc ) 


(17) 


Here, a> pe 2 = 47cn e e 2 /m e is the square of the electron plasma frequency with an electron mass m* and 
charge e, A.De 2= T e /(47m e e 2 ) is the square of the Debye length. L e = In (A) is the Coulomb logarithm. It 
has a typical value around 10 to 20. The plasma-neutral collision frequency is v en = n„ < a V t h e > . For 
typical conditions of a Hall thruster, the effect of Coulomb collision (v e i) may be smaller or of the same 
order in comparison with the plasma-neutral collision (v eil ). 13 ’ 40 


fb) Plasma- wall interactions: The inelastic electron collision with the wall allows the electrons to move 
across the magnetic field toward the anode, giving rise to “near wall conductivity’. Thus, for modeling 
the near wall conductivity, one needs to specify the secondary emission and sheath potential The wall 
with high secondary electron emission 8 can give rise to high cross-field conductivity, since a large 
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fraction of the incident energetic electrons are returned to the plasma as cold electrons with new guiding 
center drift along the direction of the electric field. 


The modeling of the sputtering yield Y and isolating its effect on the performance of the thruster is 
complicated by the fact that as the plasma energy varies along the channel, close to the inlet, wall may 
not sputter at all whereas, near the exit, sputtering yield may be considerable implying that only a 
fraction of the accelerated ions across the channel strike the wall. Based on the experimental 
observations, we shall use an empirical formula used for sputter yield 42 "* 3 , 


Y = 




H, 


(18) 


where S = lxl 0* 2 is the sputtering yield factor 37 , H s = 3000 K is the sublimation energy of boron nitride 
and Ei is the incident ion energy on the target. In the present work, we shall assume £, = 0.1 T e . For a 
near wall sheath potential, 


<P =- 


0.5 + In 


(i-*) 


( 1 - 7 ) 


m, 


y/2 


2 nm. 




(19) 


electron-wall collision frequency, for a channel of width h can be given as, 




2V„ 


'^-e v \<p£ 0 , 


2K 


( 20 ) 


the 




Here q> =e q>/T e and coefficient of secondary emission for Boron nitride wall is given as, 

1 0.576 


£ = 0.198x r 


( 21 ) 


Owing to small inertia, electron response time is much faster than the ion response time. As a result, 
electron will attain steady state much faster than ions. Keeping this in mind, electron momentum and 
energy equations are solved as steady state equations, whereas for ions and neutrals, a set of time 
independent continuity and momentum equations are simultaneously solved. The axisymmetric 
cylindrical thruster plasma is modeled by ID geometry where z corresponds to axial direction and 9 is 
along the azimuth. Following one-dimensional equations are solved in the present work. 43 
Electron momentum equation: 


oz m e n e oz m e 
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co: 


\V',+v en +a B co cJ 


V -v (V -V ) -v (V -V ) 

ez v ei\ r ez r iz J en^ez nz / 


(^-FJ+vX 


w 


where m e is the electron mass, n e is the electron number density. V ez , V a , are respective electron, 
ion and neutral axial velocities. Fe= E z /B r is the azimuthal electron drift velocity, p e = n e T e is the 
electron pressure with T e as electron temperature in eV, E z is the axial electric field, co c = eB/m e is the 
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electron-cyclotron frequency, and the source term due to ionization, recombination and charge exchange 

is S = S recom b+ Sioniz + S C ex. Following relation between azimuthal and axial velocities is utilized, 

/ \ 




a h 


V ei +V er, + (X B 0) C J 


V„ = QK„ 


(23) 


where, as is the Bohm diffusion coefficient and Q is the Hall parameter. Typical value of Hall 
parameter varies between 100 - 1000. 


We note that the suppression of axial electron mobility is due to the imposed radial magnetic field. The 
ion mobility remains unaffected by such a field. This allows plasma to support an electric field with a 
potential difference close to the applied voltage. Thus, we shall use equation (22) to determine the 
plasma potential inside the thruster. The dynamics of the electron is determined by the pressure gradient, 
by the electric and magnetic forces and by collisional exchange of momentum in equation (22). In the 
regions of sharp flow gradients, the effect of convective term may become finite and therefore, the 
convective term is retained in this formulation. Similarly, since collision time scales are much larger 
than the electron-cyclotron gyration time scale, one may ignore elastic and inelastic collision terms in 
comparison with the Lorentz force term VxB in the momentum equation. Such an approach, will 
exclude the dynamics of momentum exchange as well as the effect of ionization and recombination, 
severely limiting the applicability of the model to the thruster plasma. Furthermore, in addition to the 
presence of electron-ion and electron-neutral collisions, electron-wall collision is thought to play an 
important role in the electron transport. 13 


It is known that the classical short-range, binary collision between plasma particles v e j and plasma- 
neutrals v en are not sufficient to explain the cross field transport of the electrons and either by invoking 
Bohm diffusion 19 or by invoking plasma side-wall interaction, 13 ’ 40 such a behavior is explained. We 
model plasma wall interaction by introducing electron-wall collision frequency v w . Further, the effect of 
anomalous Bohm conductivity have been included qualitatively by including the equivalent frequency 
vg = as cd c , that incorporates the effect of magnetic field fluctuations. 


Neglecting the effect of radiation, viscous dissipation and thermal conduction, electron energy equation 
can be written as 43 


d_ 

dz 


nV 


L 


kq+tfy? s r 
I 2 2 e 


-”' eV a ^r = l—n e V' i (T i -T c )+3^ n c v m (T„ - T t ) 


dz 


m, 





(24) 


Here T e , T t and T„ (~.3 eV) are electron, ion and neutral temperatures in eV, respectively, and Ej is the 
ionization energy of the Xenon. Equation (24) includes the effect of Joule heating, contribution due to 
the exchange of random thermal energy and due to the ionization and recombination and interaction of 
the plasma with the wall. The convective flux of kinetic energy includes the flux of azimuthal electron 
kinetic energy V 2 = V ez 2 + V^ 2 . The value of a is between (2 -3). 19 
Ion continuity equation is 43 
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( 25 ) 


8n , . d ( n i V J _ o .. _ 

ir +_ 5 — 5 ~ VJh - 

In ion momentum, the momentum exchange due to collision with electrons will not be significant as ion 
mean free path is generally larger (~0.3 m) than the size of the thruster (-.02 m). Also, we consider ions 
as unmagnetized, since the gyration radius of ions is typically large for a 200G field with an ion velocity 
4xl0 3 m/s. Thus we ignore the effect of magnetic field on the ion transport. The pressure term in ion 
momentum equation can be ignored as the thermal energy of the ion is much smaller than their kinetic 
energy i.e. Then ion momentum becomes, 43 


at B & 








jV 




( 26 ) 

vU 


Neutral continuity: ^ + ^s V **l = -V . 

3 at dz 


(27) 


Equations (22)-(27) are supplemented with the current and mass conservation equations respectively as 

en i (K-F a ) = J r (28) 


mnV + m,ny. = — 

n n nz i / tz * 

A 


(29) 


Here Jt = Id/ A is the total current density; Id is the total discharge current, A is the cross section of the 
thruster channel and m is the mass flow rate. 


Before numerically solving above set of basic equations, the physical variables are normalized using 
experimental data. The mass flow rate of the propellant is th = pVA . Then the flux of the propellant is T 
= 10 23 m' 2 s' 1 . Temperature T e is normalized to first ionization potential of Xenon, T* = £,■ (12.1 eV). 
Then all dependent variable can be normalized from V* = v(7 Vm,) m/s, n*- T*!V* m' 3 , v* = a* T* s' 1 
where a* = cr 0 V(m,/m e ), oo= 3.6x 10' 20 m 2 for Xe. The fundamental length scale can be defined in terms 
of characteristic velocity and collisional frequency as, /o=F*/w The time scale is to = w' 1 . 


III. Improved Numerical Modeling 

Finite element (FE1 techniques are well known for their adaptability to arbitrary multidimensional 
geometries and accurate imposition of complicated boundary conditions. 4 4 ^ 7 Here, a multi-dimensional 
finite element formulation will be employed to solve charge-neutral continuity, momentum and potential 
equations. These equations may be expressed as L(U)=0, where U contains all state variables, e.g., ion 
density, ion velocity, potential, and L is a differential operator. The weak statement associated with a 
variational integral underlines the development of this numerical algorithm. The physical domain is 
spatially semi-discretized ( h approximated) using generic computational domain, i.e., the finite element. 
The state variables are interpolated inside the element, via the trial space FE basis set Ntfxj) that 
typically contains Chebyshev, Lagrange or Hermite interpolation polynomials complete to degree k and 


13 



to dimension j. The spatially semi-discrete FE implementation of the weak statement WS* for L(U)=0 
leads to 


WS" = S. 






= 0 


(30a) 




\N^\N t (s)dT + |J W,(f r f;) V<T (30b) 

va n e n, an, nan* J 

5 e symbolizes the “assembly operator” carrying local (element e) matrix coefficients into the global 
arrays, s is a source term (e.g. inelastic processes) and f} and if are convective and dissipative flux 
vectors, respectively, and h is the direction normal. Application of Green-Gauss divergence theorem 
“weakens” the order of derivatives in (30a) and yields natural homogenous Neumann honndarv 
conditions. The surface integral in the second line (30bl contains the (uni known boundary fluxes 
wherever fixed or flux boundary conditions are enforced accurately. 


•dN„ 


Independent of the physical dimension of the working domain G, and for general forms of the flux 
vectors, the semi-discretized weak statement always yields an ordinary differential equation (ODE) 
system that is fully discretised using a ^-implicit or r-step Runge-Kutta type time integration. The 
terminal ODE is usually solved using a Newton-Raphson scheme for U(t): 48 " 49 

U™ =u; + , +AIT =U, +y.U p+1 , where AU' = - [M + 6M(8R / 5U)]" 1 R(U) (31) 

p = o 

Here, a ^-implicit time marching procedure is employed. In (31), M = 5 e (Me) is the “mass” matrix 
associated with element level interpolation, R carries the element convection, diffusion and source 
information. The calculation of the “Jacobian” 8Rld\J and inversion of the M+6Ar(5R/5U) matrix with 
sufficient accuracy is obviously a numerical challenge. However, unlike the traditional finite 
difference/volume methods, the present FF. algorithm allows one to simulate the system simultaneously 
without requiring any sub-iteration for the Poisson solver. 

The stiffness of the jacobian matrix for such coupled nonlinear problem becomes very high for realistic 
mass ratio of electron and ion. This results in solution divergence for standard finite 
difference/volume/element numerical approach on a moderate to fine mesh and require filtering and 
other expensive approaches for solution control. As a remedy, we utilized an efficient high-order 
accurate, sub-grid embedded (SGM1 finite element method 50 ' 51 to achieve stable monotone solution on a 
relatively coarse mesh (see Figure 2). SGM elements utilize a non-linear sub-grid basis function, which 
is determined based on the local fluid velocity, mesh size and conductivity. SGM also serves as an 
excellent preconditioner for parallel solvers as it improves positive definiteness in the element matrix. 
Further details of the code are described in Ref. 44-46 and in Ref. 49. 
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Figure 2. Coarse mesh comparison between standard Petrov-Galerkin and sub-grid embedded solutions 
for the ion number density distribution after 10" 6 sec (from Ref. 44). 


The choice of time step is dictated by the Courant-Fredrich-Lehy condition. 52 The solution at any time 
step is declared convergent when the maximum residual for each of the state variables becomes smaller 
than a chosen convergence criterion of e=10 4 . The steady state is declared when the above 
convergence criteria is met at the first iteration of any time step. Here, the convergence of a solution 
vector U on node j is defined as the norm: 


EizM< e 

IlUyll 


(32) 


IV. Numerical Results and Discussion 

We recall that the thruster plasma is modeled by a 2D axisymmetric (r-z) geometry, where r corresponds 
to the radial direction and z corresponds to the axial direction. The 0-direction is along the azimuth. We 
consider a two-dimensional magnetic field with radial and axial components, where the radial field is 
dominant (Fig. 3). The magnetic field lines near the exit close outside the thruster indicating that the 
near plume region plasma will be affected by the presence of the magnetic field. However, the 
simulation domain in the present work corresponds only to the bounded region. We have employed a 
trapezoidal time stepping procedure, i.e. >9=0.5 in (31), for this model. In the present formulation the ion 
dynamics are time dependent, whereas electron dynamics have been assumed time-independent. This is 
a plausible assumption since, owing to their small inertia, the electron will reach a steady state over the 
ion dynamic scale. 10 The code uses variable time steps until the transient features die down and the 
iteration converges to a steady state. All solutions presented in this section have iterated to a steady 
state. 
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Figure 3. Measured radial magnetic field in Gauss inside the thruster acceleration channel. 


Fig. 4a describes the neutral density contours. We use the reference values of physical parameters 
pertaining to a typical 1.5 kW class thruster that has a mass flow rate ~ 6 mg s' 1 corresponding to n„ = 
10 18 m' 3 . As is clear from the figure, neutral density is the highest at the inlet region and gradually 
decreases towards the channel exit. As ionization increases towards the exit (due to an increase in the 
electron temperature'), the neutral number density should decrease and we see such a behavior. This is 

A / / * 

consistent with the fact that as a neutral enters the thruster chamber it undergoes impact ionization. 
Some experiment in the literature 10 suggests that the minimum in neutral density is not necessarily 
correlated to the corresponding increase of the ion density. However, in the present work the plasma 
density prediction displays a correlation between ion and neutral number densities that is similar to the 
reported experimental data. 2 ' 6 We attribute this correlation to the temperature dependent, self-consistent 
calculation of the ionization rate. 


Fig. 4b plots the plasma number density contours. The ion (electron) number density increases rapidly 
from a base value of 10 17 m' 3 and attains a maximum value 7xl0 17 m' 3 upstream of the acceleration 
channel before decreasing near the exit. The experimental results 2 ’ 14 show that the plasma density 
reaches its peak value inside the acceleration channel, near the inner wall before the exit plane. In this 
region, the radial magnetic field is the maximum and thus a large number of electrons are inhibited from 
moving in the axial direction, resulting in a high probability of impact ionization and plasma production. 
The maximum plasma density inside the acceleration channel agrees with the fact that the ionization 
channel is well inside the thruster. There is no significant effect of ionization and recombination on the 
plasma number density. This could have been anticipated on the grounds that in a Hall thruster, where 
the pressure is low and ion currents exceed the electron current, the effect of the ion production and loss 
to the ion continuity equation (6) is negligible. 20 
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(a) 


Figure 4. (a) Neutral number density contours in m' 3 , (b) Plasma number density contours in m 


(b) 
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The experimental result for a 1.6 kW class thruster 13 displays two distinct peaks in the ion number 
density profile located at about 0.02 m and 0.032 m from the anode. These peaks are attributed to 
different ionization mechanisms - to the electron thermal energy upstream (0.02 m) and to the 
availability of electron gyration energy at 0.032 m. These results underline the complexity of the 
thruster plasma dynamics and inadequacies of the existing numerical models. Several important 
questions need to be addressed in order to explain the physical mechanism behind the experimentally 
observed transition from double hump to single hump ion density profile when operating at 1 .6 kW and 
3 kW. 14 If at 1.6 kW, plasma undergoes a unique ionization-recombination-ionization cycle, then such 
behavior should be reflected in the neutral velocity and density profiles. We anticipate that at higher 
operating powers of the thruster, the neutral 
velocity will exhibit an initial increase 
(corresponding to the loss of slow neutrals due to 
ionization, i.e. number of fast neutrals increase), 
then decrease and again an increase. Also, the = 0 064 E" 
neutral number density should exhibit an initial 
decrease, then an increase and again a decrease in 
its number density. It points to the necessity of 
generalizing the numerical model on the one hand 
and further experimental investigation of the 
thruster plasma dynamics on the other hand. 
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The direction of ion streamlines as plotted in 
Figure 5 shows that the ion flow diverges towards 
the sidewalls in the downstream section of the 
channel, indicating the presence of a radial electric 
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Figure 5. Directed ion trajectories show the 
decollimation of ions near the channel exit. 


field. Haas has experimentally inferred the 

presence of this radial field, where the radial asymmetry in the ion number density has been attributed to 
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the presence of such a field. Figure 6 shows that the magnitude of the radial velocity contours increases 
in the region of strong magnetic field (B r ) confirming the experimental observation. 2 The interaction of 
ions with the ambient magnetic field could be another possible reason for divergence of the ion 
streamlines. The magnetic field, which confines the electrons in the azimuthal direction and inhibits 
their axial motion, may exert its influence on ions due to the collisional coupling of ions with the 
electrons - like ambipolar diffusion in the interstellar plasma. 53 Thus, even though ions are not directly 
coupled to the magnetic field, they may interact with the magnetic field through the electrons. We also 
anticipate that an additional divergence in the ion- 
beam may appear once the thermal pressure 
gradient is included in the ion dynamics. When the 
ion pressure gradient is taken into account, it will 
give rise to a radial electric field that may cause a 
dispersion of the accelerated beam. The 
decollimation of the ion beam in the radial 
direction will reduce the thrust. 


In addition, in the presence of ionization and 
recombination there is a slight increase in the ion 
radial velocity near the exit as shown in Fig. 6. 
This increase in the radial velocity could be due to 



mv 
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of the slow ions to the walls. 


Distance from Anode (m) 

Figure 6. Radial velocity contours in m/s. 

However, a definite correlation between the 

velocity and density can only be made if plasma-wall interactions are also included in the model. 
Recent numerical studies 21 suggest that plasma-wall interaction may affect the plasma density, near the 
exit plane. 

Figure 7 describes the electron temperature contours. 

We note that the increase in die temperature is not 
uniform in the channel. The maximum increase of 
~28 eV occurs just downstream of the center of 
channel, towards the inner wall. The peak in the 
temperature can be attributed to the Ohmic heating 
due to the maximum gyration energy in this region. 

The rise in temperature is similar to the measured 
electron temperature near the exit of the 1.6-kW class 
thruster of Hass and Gallimore. 14 However, the 
experimental electron temperature peak is spread 
along a radial line concentrated near the channel exit. 

Our numerical electron temperature results do not 
clearly reproduce this profile and may point to the 
limitation of the present model. Secondary electron 

emission, ion sputter yield, the electron near wall conductivity, the near wall sheath effect etc., may all 



Figure 7. Electron energy distribution in eV. 
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affect the electron temperature profile. The present model will be extended to include plasma-wall 
interactions for better benchmarking with the experimental data. 

Figure 8 shows the potential contours inside the 
acceleration channel. The potential is highest at 
the inlet (near the anode) and is set to zero at the 
outlet (near the cathode). This is similar to the 
numerical approach presented by Haas. 2 
Although the computed potential is set to zero at 
the channel exit, observations 2,14,17 indicate that 
only one half to one third of the potential drop 
actually takes place downstream of the thruster. 

In this numerical model, the full potential drop is 
forced to occur inside the channel. We further 
note that the potential inside the channel is very 
sensitive to the boundary conditions on the 
plasma velocities. Such effects will be evaluated 
in more detail if future support is available. 

In Figure 9, the plasma number density, electron temperature, plasma potential along with the electric 

field and electron gyration energy are given at radial cross section r = 0.056 m, centrally located 

between the outer and the inner wall. We see that the plasma number density exhibits a narrow half 

width at full maximum (hwfin) just upstream of 2 = 0.1 m, suggesting a localized ionization region. It 

was noted by Kim 5 in one of the early experiments on SPT that a very small ionization region precedes 

the acceleration zone. Our numerical result is consistent with the experimental observations. 2,5 The 

electron temperature T e profile in Fig. 9 predicts a 

sharp increase near the inlet and then a gradual rise 

before reaching the maximum of about 22 eV near 

the two-third length downstream of the channel. 

This is not consistent with the trend in plasma 

density profile. One would expect that the plasma 

density maximum will coincide with the electron 

temperature maximum, and that the first local 

maximum in electron temperature at around 14 eV 

at normalized axial location of 0.11 will be 

reflected in a small peak in the plasma density. 

However, the first local maximum in the 

temperature may reflect the loss of slow electrons 

due to recombination, since plasma number 

density is not a sensitive function of source or sink 

term 20 , the local electron temperature maximum of 

14 eV (at z/L = 0.1 1) does not affect the plasma density profile. ** - 

' r density, and temperature below the centerline. 
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Figure 8. Potential distribution inside the 
acceleration channel in volts. 
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The computed potential profile is consistent 
-3V/+6V) for a 1.5 kW thruster that operated 
at 300 V. 16 The potential profile is flat in 
most of the channel and approaches zero at 
the exit. Since potential is forced to zero at 
the exit, the simulation result starts diverging 
with the experimental result near the exit. 
The gyration energy (last frame in Fig. 9), 
displays an increasing trend, reaching 
maximum near the exit, which agrees well 
with the published Hall contours 2 ’ 21 that 
display a maxima upstream of the channel 
exit for 1.6 and 3 kW-class thrusters. Finally, 
the thrust at the exit plane of the thruster 
acceleration channel is calculated via Eqn. 
(24) of Haas and Gallimore. 54 The simulation 
result shows a steady state thrust of 79.4 mN, 
which is within the measured data of 95 3 
mN and the calculated value of 68 mN at the 
exit plane of the 1 .6 kW thruster. 54 


with the known experimental data (with uncertainty of 
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Figure 9 (cont.). The distribution of potential and 
electron gyration energy below the centerline 
(r=0.056m) of the acceleration channel. 


In Figs. 10-13, the solutions of (22)-(29) are also documented to show the effect of sputter and SEE. The 
applied magnetic field has a peak upstream of the exit plane. The dotted line correspond to a case when 
only effect of secondary emission is considered in the plasma-wall interaction and bold line correspond 
to a case when both secondary emission and sputtering yield effect have been considered. The plasma 
number density (Fig. 10a) increases rapidly from a base value of 2.8xl0 17 m* 3 and reaches a maximum 
value 1.6xl0 18 m' 3 upstream of the acceleration channel before decreasing near the exit in the presence 
of secondary emission. However, when both secondary emission and sputter yield, due to ion-wall 
interaction are included, the decrease in number density is more pronounced toward the exit (bold line). 
The change in plasma density toward the exit in the presence of sputtering yield is consistent with the 
fact that plasma is partly lost to the wall and the effect is more pronounced near the exit of the channel. 
The experimental results 19 shows that the plasma density reaches its peak value inside the acceleration 
channel, right bottom comer of the exit plane. The location of ionization zone (~0.6) is same for both the 
cases. The maximum plasma density inside the acceleration channel is in agreement with the fact that 
the ionization channel is well inside the thruster. 


The rapid increase in the ion number density is reflected in the rapid decrease in the neutral number 
density (Fig. 10b) from 2xl0 18 m' 3 to approximately 1.6xl0 18 m' 3 . This is consistent with the fact that as 
neutral enters the thruster chamber it undergoes the impact ionization. The effect of sputter yield and 
secondary emission (bold line, Fig. 10) is not very significant in comparison with the case when only 
secondary emission is present (dotted line). 
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Figure 10. Effect of SEE and sputter yield on (a) ion and (b) neutral number densities. 

Fig. 11a describes the axial ion velocity profile. The velocity peaks downstream of the channel, before 
the exit. This indicates that the location of the acceleration channel is inside the acceleration channel at 
0.75. The ion velocity is slightly higher in the presence of sputtering (bold line) than in the presence of 
secondary emission only (dotted line) though difference is not very large. Ions are accelerated mainly 
due to the presence of the potential gradient, which is maximum near the channel exit, Fig. 6. Further, 
one may infer from the location of the acceleration channel that die width of the ionization region is 
narrower (-0.15) than the width of the acceleration channel (-0.25). This is in conformity with the 
experimental results. 6 ’ 25 

Fig. lib shows the electron velocity profile. There is a slight increase in die electron velocity due to the 
slight increase in the electric field (Fig. 12, Bold line) in the presence of sputtering effect on the plasma- 
wall collision frequency. The electrons from the cathode, located just outside the chamber of a Hall 
thruster, are accelerated towards the anode. Large negative velocity near the exit is consistent with the 
large electric field, which are responsible for accelerating the electrons towards the inlet. These inward 
moving electrons, on their way to anode collide with the neutrals and ionize them. As a result, electron 
velocity decreases towards the anode as reflected in the figure. 



Figure 1 1. The effect of HET wall interactions on axial (a) ion and (b) electron velocities. 

Figure 12 shows the potential profile inside the acceleration channel. The change in potential profile is 
not significant in two cases. We see that the potential has a zero gradient inside the thruster channel 
similar to the experimental data. 6 However, the computed potential vanishes at the channel exit, while 
observations 6,17 indicate that only one half to one third of the potential drop takes place downstream of 
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the thruster exit. This difference is due to the imposition of zero potential boundary condition at the exit 
plane in numerical simulation, i.e., full potential drop is forced to occur inside the channel. 



Figure 12. Electric field E and potential difference 4> - <)>e. The potential remains unchanged for 2/3" 1 of 
the channel and then sharply drops to the exit potential <t> E . 

In the presence of sputtering effect, azimuthal electron drift velocity (Fig. 13) profile towards the exit is 
quite different than when only secondary emission is present in the plasma-wall collision frequency. The 
azimuthal velocity V e g increase in the presence of sputtering yield is consistent with the electric field 
profile (Fig. 12 bold line). Towards the exit, V e6 is smaller than in the absence of sputtering (dotted 
line). This behavior indicates that the plasma- wall interaction affects the potential towards the exit. 
This is consistent with the change in electric field profile in file previous figure. The drift velocity is a 
consequence of the crossed electric and magnetic field and gives rise to Hall current density, JH»en e Vo 
. In Fig. 13, the peak in the azimuthal velocity considering sputter yield (bold line) is consistent with the 
electron temperature profile. 



Figure 13. Electron drift velocity with sputter yield predicts the maximum just upstream of the channel 
exit (bold); the trend is different in the case without the sputter (dotted line) where the peak is 
at the exit plane. 

The temperature distributions in two cases are not very different. The increase in the temperature is not 
uniform in the channel. The maximum increase occurs just downstream of the center of the channel in 
both cases. However, when plasma-wall interaction considers the ion sputter yield effect, electron 
temperature decreases slightly towards the exit (bold line). The peak in electron temperature can be 
attributed to the Ohmic heating due to the maximum gyration energy in this region. This trend in 
temperature distribution is similar to the results reported in the literature. 19 The computed temperature 
profile is in agreement with the measured electron temperature near the exit. 6,25 
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Finally, the new NASA-457M (Fig. 14) is modeled and simulated for two flowrates and for two 
different fluids, namely, xenon and krypton at 464 seem and 829 seem. The 464 seem case corresponds 
to an anode mass flow rate of 29.1 mg/s for krypton and 45.8 mg/s for xenon. The 829 seem case 
corresponds to an anode mass flow rate of 52 mg/s for krypton and 86.4 mg/s for xenon. For a given 
proprietary magnetic field and a fixed channel (2.5"x2"), four cases are simulated by imposing a 
potential drop across the channel. Table 1 summarizes these results. 



(a) Xenon glow (b) Krypton glow 

Figure 14. NASA-457M during testing. 


Gas 

Discharge 

Voltage 

(Volts) 

Anode 
Massflow 
rate (mg/s) 

Calculated 

Thrust 

(mN) 

Specific 

Impulse 

(seconds) 

Discharge 

Efficiency 

Xenon 

400 

45.8 

989.4 

2201 

59.2% 

86.4 

1983 

2240 

61.3% 

Krypton 

400 

29.1 

703.2 

2463 

53.2% 

52.0 

1335.1 

2617 

56.3% 


Table 1. NASA-457M computed results. 


V. Conclusions: 

In this report, a finite element based multidimensional formulation of partially ionized plasma using 
multi-component fluid equation is developed in the presence of ionization, recombination and wall 
losses. The effect of inelastic processes viz. secondary emission and sputter yield has been incorporated. 
The model is then applied to study the dynamics inside the Hall thruster. The ion and neutral dynamics 
are time dependent, whereas electron dynamics is assumed time-independent. By using a third order 
electron temperature-dependent polynomial, a self-consistent calculation of the ionization rate has been 
carried out in the model. Our simulation results suggest that the increase in the plasma number density 
is correlated with the decrease in the neutral number density. The plasma density prediction is similar to 
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the reported experimental data. 2,5 There is no significant effect of the ionization and recombination on 
the plasma number density. This fact is consistent with reported observation 20 . 

The electron temperature inside the channel shows a gradual increase at the centerline, and predicts a 
hump upstream of the exit at a location between the centerline and the inner wall, which is in agreement 
with the experimental observation 2 that shows a peak next to the exit for a 1 .5 kW class thruster. The 
potential profile agrees with recent experimental studies. 17 The ion streamlines suggests that ions are 
primarily accelerated due to the axial electric field and reach the maximum velocity near the exit plane. 
The ions are moving towards the sidewalls near the thruster exit, indicating the presence of a finite 
radial electric field. Experimental data 2 confirms the presence of such a radial electric field. The 
numerical results are very sensitive to boundary conditions. Therefore, boundary condition issues 
require a detailed investigation. Also, proper modeling of the plume region is necessary to make this 
model useful for design purposes. 
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